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Abstract 


This work presents a numerical model for aggradation and degradation in tree-type 
channel networks. The model proposed numerically solves the differential equations 
for water flow and sediment flow taking the quasi-steady flow assumption and using 
uncoupled technique. The model is used to simulate aggradation and degradation 
profiles in converging and diverging tree-type networks. Both rectangular and non- 
rectangular cross-sections are considered. Three different procedures for changing 
the cross-sectional shape of a channel are adopted. It also includes (1) Power law 
model, (2) Iowa model and, (3) Meyer-Peter-Muller equations for computing the 
sediment transport capacity. Applicability of the model is demonstrated through 
several simulation studies for (1) aggradation due to sediment overloading and, 
(2) degradation due to sediment depletion. The numerical results shows that the 
proposed model satisfactorily simulates the basic features of aggradation and degra- 
dation process. 
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Chapter 1 


Introduction 


1.1 Overview 

Most natural beds are composed of loose materials which can move due to the 
bed shear stress exerted by the flowing water. The natural long-term (geological 
time scales) evolution of these channels known as alluvial channels, produces slopes, 
widths, depths and velocities such that the flows are able to transport the imposed 
sediment discharges. The amount of sediment coming into a particular reach of 
the channel is equal to the sediment going out from that reach and therefore, the 
stream bed elevation will not change over a long period of time. Such alluvial 
streams are said to be in regime or in equilibrium. Any short term, of order of hours 
to weeks, changes in flow and sediment transport rates from long-term values are 
accommodated by the channel by adjusting its depth, velocity and the roughness. 
On an intermediate time scale, of the order of months to years, channels can also 
adjust their large-scale geometry i.e, widths, depths, bed levels and plan form, in 
response to natural and, or human interference in their hydrologic and sedimentary 
regimes. 

Many times, alluvial channels are harnessed for power generation, navigation, 
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and the interdependencies among them in river mechanics are too complex to permit 
simple analytical solutions. River response to the structural modifications extends 
over long distances and over many years. Therefore, evolution of this response 
using physical model studies is not feasible. The only alternative left is the use of 
mathematical models, however rudimentary they may be. The objective of the study 
reported herein is to develop a numerical model for aggradation and degradation of 
tree-type channel systems whose equillibria have been disturbed. Literature in the 
are of movable bed modeling is reviewed in the following section to place the present 
work in perspective. 

1.2 Literature Review 

Analytical solutions have been developed of simplified governing equations describ- 
ing complex phenomenon of agradation and degradation processes. De Vries (1973) 
suggested that the flow can be considered quasi steady for the study of bed level 
variation. A linear parabolic equation for the bed level variation was obtained by 
assuming the flow to be uniform, and steady in a wide rectangular channel. Soni 
et al. (1980) used this model to predict the transient bed profiles due to sediment 
overloading. Jain (1981) pointed out an error in their boundary conditions and 
presented an analytical solution utilizing more appropriate boundary conditions. His 
computed results compared satisfactorily with the experimental data. Begin et al. 
(1981) used a similar model to compute longitudinal profiles produced by base levels 
lowering. Gill (1983a, 1983b) solved the linear diffusion equation for agradation and 
degradation by Fourier series and by the error function methods. Jaramillo and 
Jain (1984) developed a non-linear parabolic partial differential equation and solved 
it by the method of weighted residuals. They compared their computed results 
with experimental data obtained by Newton (1951) and Soni et al. (1980). Zhang 
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and Kahawita (1987) detected an error in the solution of Jaramillo and Jain (1984) 
and presented a numerical solution of the non-linear parabolic model. Gill (1987) 
also presented non-linear solutions for aggradation and degradation . Ribberink et 
al. (1987) suggested that a non-linear hyperbolic model should be preferred to a 
non-linear parabolic model. 

Analytical models for aggradation and degradation are based on the assumption 
of steady water flow. This assumption is not valid while estimating the bed level 
changes during unsteady water flow conditions. It is also strictly not valid even if 
the discharge is constant since the water surface profile computations depend upon 
the bed slope. Analytical models assume also a simplified equation for sediment 
transport capacity and a wide rectangular cross section. Therefore, water flow 
equations along with sediment continuity equation are generally solved by numerical 
techniques. Since the early work of Vreugdenhil and De vries in 1967 (Cunge et al. 
1980), several one dimensional movable bed models have been developed. Holly 
(1986), Dawdy and Vannoni (1986), Cunge et al. (1980) and Palaniappan (1991) 
presented excellent reviews of numerical simulation of alluvial hydraulics. Lu and 
Shen (1986) tested several numerical aggradation and degradation models by com- 
paring the computed results with the laboratory data obtained by Suryanarayana 
(1969). Most of the standard one dimensional unsteady sediment transport models 
can be classified into either uncoupled models, wherein the water flow equations 
and sediment continuity equation are uncoupled during a given time step or coupled 
models, wherein all the governing equations are solved simultaneously. They can 
also be classified based on whether full Saint Venant equations for shallow water 
flow are considered (Unsteady flow models), or the water flow is assumed steady 
during the computations of bed level variation (Quasi steady flow models). 
Quasi-steady flow uncoupled models : HEC-6, developed by US Army corps 
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of Engineers (Thomas and Prashun 1977, Thomas 1982, Copeland 1986), Flvial- 
3 (Chang and Hill 1976), KUWASER (Simons et al. 1979), IALLUVIAL (Karim 
and Kennedy 1982, Holy et al. 1984, Karim 1985, Holly and Karim 1986) and a 
recent model developed by Palaniappan (1991) are good examples of quasi-steady 
uncoupled models. Witkowska (1971), Bettess and White (1979, 1681), and Young 
(1984) also developed similar models. In all these models, one- dimensional steady 
flow equations are used to compute the gradually varied flow profile in a channel. 
The sediment mass balance equation is then used to determine the bed level change 
during the specified computational time step. This is followed by recomputation 
of the water surface profile corresponding to new bed levels, and the procedure 
is repeated. Models belonging to this class differ from one another in terms of 
(1) the equation for the sediment transport capacity, (2) the equation used for 
computing the friction slope, (3) the finite difference approximation of the sedi- 
ment mass balance equation, (4) the level of detailing the cross-sectional shape 
and, (4) the armouring procedure. For example, Witkowska (1971) ignores the 
armouring procedure completly while Holly and Karim (1986) include, probably the 
most sophisticated procedure for simulating the armouring phenomenon. Recently, 
Mohapatra and Bhallamudi (1994) used a quasi-steady flow model for studying the 
be level variation in channel transitions. 

Quasi-steady Flow Coupled Models : CHAR-3 (Cunge and Simons (1975) and 
CHAR-2 (1978) belong to this class of mobile bed modeling. The steady flow 
momentum equation for water and the unsteady mass balance equation for the 
sediment transport equation are solved using the Preismann implicit scheme. The 
resulting nonlinear algebraic equations are solved simultaneously to achieve coupling 
between the water flow and the sediment movement. CHAR-2 and CHAR-3 are 
based on basic modeling philosophy proposed by Cunge and Perdreau (1973). 
Unsteady Flow Uncoupled Models : Quasi-steady flow models can not be 
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applied strictaly while studying flood related changes. Complete Saint Venant 
equations for water flow should be considered in such cases. FLUVIAL programs 
developed by Chang (Chang 1982, 1984, 1985), the model developed by Bennet Nordin 
at USGS (1977), the models developed at the Colorado State University (Chen 1973, 
Chen et al. 1975, Chen and Simons 1980) and CHAR-1 (Cunge et al. 1980) are 
examples of unsteady, uncoupled models. The water flow and the sediment transport 
computations lag by one computational time step in all these models. 

Unsteady Flow Coupled Models: Vreugdenhil (1982) and Lyn (1987) studied 
the stability characteristics of one dimensional movable bed modeling. Lyn (1987) 
used perturbation techniques to identify multiple time scales in the governing equa- 
tions and suggested that complete coupling between the full unsteady water flow 
equations and the sediment continuity equation is desirable in cases where the 
conditions are rapidly changing at the boundaries. The models developed by Chen 
and Simons (1975), Chollet (1977) and Krishnappan (1981,1985) belong to the class 
of one dimensional, unsteady, coupled models. Park and Jain (1986) achieved a 
partial coupling between sediment flow and water flow by iterating the solution 
whenever the spatial gradient of change in bed level became too large in order to 
avoid numerical instabilities. Later, Park and Jain (1987) introduced a simple bed 
layer model to simulate the armouring process. Lai and Chang (1987) developed an 
unsteady model using characteristic method with implicit and reach back schemes. 
Havno et al. (1985) modified the model MIKE-11 developed at the Danish Hydraulic 
institute (DHI 1983) to include the erosion procedures. In this model, the St.Venant 
equations are solved using a fully centered implicit difference scheme. Rahuel (1988) 
and Rahuel et al. (1989) developed a fully coupled implicit finite difference model 
which includes a procedure for routing the graded mixtures. This model requires 
the solution of (2k+l) partial differential equations, where k is the total number of 
sediment size classes. Later, Holly and Rahuel (1990) extended this model to treat 
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suspended and bed load transports separately. It should be noted that this model 
was developed for a single channel only. It requires enormous data preparation 
and computational power for application to even single channels. Extensions of 
this model to channel networks is a very difficult task. Correia et al. (1992) also 
developed a fully coupled unsteady mobile boundary flow model for application to 
a single channel. In this model, the governing equations are formulated in such 
a way that water flow equations and sediment transport equation are explicitly 
coupled. The equations are solved using the Priessman scheme. However, a new 
method is introduced for solving the simultaneous linear equations which result 
from the discritization. It should be noted here that most of the unsteady, coupled 
models use implicit schemes to solve the governing equations. Although, large 
computational time steps can be used with these schemes, they involve solution of a 
system of equation by matrix inversion during each computational time step. This 
results in high computational costs, especially when extended to the case of channel 
systems. Complete unsteady flow modeling is usually required only when studying 
the flood related changes. Computations in these cases may have to be performed 
for time periods of at the most for two to three months. Therefore, use of explicit 
methods may not be a bad choice. Bhallamudi and Chaudhry (1991) presented 
one such completely coupled unsteady flow model using the second order accurate 
Maccormack scheme. This model has been used to study (1) the aggradation due 
to sediment overloading, (2) the degradation due to base level lowering and, (3) the 
knick-point migration in laboratory channels. 

Although the present trend is to use ’coupled models’ for reverine studies, some 
researchers feel that the concerns in the literature about the use of coupled models 
may have been overstated. Recently, Cui et al. (1996) compared the numerical 
results obtained using coupled and quasi-steady uncoupled models for cases with 
Froude number very close to unity and also for cases in which the upstream water 
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discharge, sediment feed rate and downstream water elevation varied strongly. There 
was a very good agreement between the numerical results obtained using the two 
models, although uncoupled models are inherently more unstable than the coupled 
models. Performance of uncoupled models appears to be satisfactory for most cases 
of riverine problems. It should be noted that quasi-steady flow uncoupled models 
require significantly less computational power for their application. Therefore, they 
are ideally suited for application to large river network problems and for simulating 
long term effects. It is also easy to include sophisticated relations for sediment 
transport capacity, armouing process and friction factor in these models. Therefore, 
a quasi-steady flow uncoupled approach is followed in the present study. 

Channel Network Models : There has not been much research in the area 
of sediment routing in channel networks. In fact, studies on even a basic topic 
such as the gradually varied flow (GVF) computation in channel networks are 
few and recent. Shulte and Chaudhry (1986) presented a numerical algorithm 
for computing water-surface profiles in steady-state, gradually varied flows in open 
channel networks. In this procedure, channel in the network is discritized into finite 
difference nodes and the energy equation is written for each pair of two consecutive 
finite difference nodes. These equations along with the compatibility equations 
for the channel junctions are simultaneously solved using the Newton-Raphson 
technique. This method can be used for both tree type and looped networks. 
However, the structure of network is not exploited in this algorithm to achieve 
computational efficiency. Also, the convergence of the Newton-Raphson technique 
is not guaranteed. 

Many researchers suggest the use of false transient approach to obtain steady 
state solutions in the networks (Holly et al. 1985). Research interest in the area of 
unsteady flow modeling in channel networks has been revived in recent years. The 
commercial model FLDWAV (Fread 1985) is applicable for river networks. However, 
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it does not exploit the special structure of tree type networks and therefore, is 
computationally very intensive. Application of implicit methods for unsteady flow 
in channel networks becomes very difficult because the double sweep algorithm can 
not be applied in a straight forward manner. Akan and Yen (1981) applied an 
iterative technique in routing floods through tree type channel networks. However, 
they used the governing equations for diffusion routing and therefore, the backwater 
effect of a junction can not be considered. If the down stream backwater effect is 
significant, a large number of iterations may be necessary to achieve the desired 
accuracy. Tucchi (1978) and Barkau et al. (1989) applied the skyline algorithm 
for computing unsteady flow in channel networks. The skyline algorithm requires 
a large number of iterations during the reduction pass. In recent years, there has 
been an emphasis on solving the entire network as a unit system rather than as 
a combination of independent systems. Choi and Molinas (1973) introduced one 
such direct solution algorithm based on the double sweep method. This method is 
applicable to dendritic channel networks. The model proposed by Choi and Mohinas 
(1993) requires operation on a 2N x 4 matrix as compared to the operation on a 
2N x 2N matrix (N = number of cross sections used) in the direct application of 
the Preismann scheme. In this model, the channel network problem is converted 
into a single channel problem by using separate recursive equations for each of the 
components of the network i.e, (1) interior channel cross sections, (2) converging 
channel junction and, (3) diverging channel junction. Nguyen and Kawano (1995) 
proposed an alternative algorithm in which all the channel sections and junctions 
can be uniformly treated. The model proposed by Nguyen and Kawano (1995) is 
applicable to a broader range of dendritic type of networks. Recently, Tiwari (1996) 
slightly modified the algorithm presented by Nguyen and Kawano (1995) to take 
into consideration Log Pearson Type-Ill hydrograph as the inflow boundary and 
a barrage as the downstream boundary conditions. The algorithm was successfully 
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applied to a channel network as large as the Ganga system from Hardwar to Farakka 
with all its tributaries. The algorithm was applied to simulate flows for one month 
during the flood season. It should be noted here that algorithms discussed so far are 
applicable for simulating unsteady flow in dendritic systems only. On the other hand, 
Kutija (1995) applied graph theoretic concepts and computer algebraic procedures to 
develop an algorithm for computing unsteady flow in combined dendrite and looped 
channel networks. The complex algorithmic structure proposed by Kutija (1995) is 
required only for application to large urban drainage systems where channels are 
multiply connected. Most of the natural river systems are tree type except in the 
deltaic region where a looped structure may exist. Therefore, algorithms which 
exploit this tree structure are most useful in the present context. It should also be 
noted here that all the algorithms for unsteady flow in networks discussed so far 
are based on the implicit finite-difference method. On the other hand, Misra et. 
al. (1989,1992) developed an explicit numerical model, TRANS-NET for study 
of unsteady flow in irrigation canal networks. Space line interpolation explicit 
characteristic scheme on rectangular grids is used for solving the Saint Venant 
equations. This model can be also used for flood routing in river networks. 

All the network models described so far are applicable only for flood routing. 
Network models for sediment routing are very few. A commercial model called 
UUWSR was developed at the University of Colorado ( Chen and Simons, 1980 ) 
for routing floods and sediment through the channel networks. This model is based 
on the complete unsteady flow equations but follows a decoupled approach. Holly et. 
al. (1985) and Yang (1986) also developed a network model for flood and sediment 
routing in multiply connected channels. Their model is based on the quasi-steady 
assumption for the water flow and uses the uncoupled approach. The quasi-steady 
flow conditions, however, are obtained using the unsteady flow equations in a false 
transient approach. It should be noted here that use of a steady flow network solver 
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in a quasi-steady flow uncoupled sediment routing model would result in significant 
saving of computational time. The steady flow network solver, recently developed 
by Naidu et. al. (1997) is ideally suited for this purpose. Recognising that the tree 
type networks are more common that the looped networks, Naidu et. al. (1997) 
developed a very efficient algorithm for computing GVF profiles in tree type channel 
networks. This algorithm is based on some basic graph theoretic concepts and the 
shooting method for the boundary value problems. It was shown that the proposed 
method could result in significant savings in computational costs, some time as high 
as 200 times when compared with the steady-flow solver developed by Shulte and 
Chaudhry (1986). It should be noted here that the algorithm developed by Naidu 
et. al. (1997) is designed for computing GVF profiles in diverging channel networks 
where water from a main channel gets distributed to branch channels. Therefore, 
a sediment routing model based on the network solver developed by Naidu et. al. 
(1997) would be ideally suited for evaluating aggradation and degradation in (1) 
channels in deltaic region and, (2) unlined irrigation canal systems. This method is 
not efficient for converging channel systems where the flows from many tributaries 
join together to form the main river. However, computation of GVF profiles in 
converging channel systems is not complicated since it is possible to determine the 
flow in each channel as the solution to an initial value problem. 

1.3 Objectives of the Present Study 

The literature review presented in the previous section shows that the majority of 
the aggradation and degradation models are applicable for single channels. Only few 
models have been developed to study the bed level variations in multiple channel 
systems. Objective of the present study is to develop quasi-steady flow, uncoupled 
models for aggradation and degradation in tree-type channel networks. Two separate 
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numerical models are developed for the cases of converging and diverging channel 
systems, respectively. Very efficient numerical procedures are adopted while develop- 
ing these numerical models. Also, general cross-sectional shapes and a sophisticated 
predictor for the sediment transport capacity are considered. This work is part a on 
going effort at 1. 1. T. Kanpur (supported by the Ministry of Water Resources, India 
) to develop a comprehensive mathematical modeling system for flood and sediment 
routing in channel networks 

The following chapters present the theoretical considerations and, numerical 
procedures . Applicability of the model is demonstrated in chapter 4. The thesis 
concludes with a summary and recommendations for further work. 
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Chapter 2 


Theoretical Considerations 


Simulation of aggradation and degradation in open-channels with movable beds 
involves the numerical solution of governing differential equations for water and 
sediment flow. Inspite of extensive research to understand the relationship between 
water flow and sediment transport (Shen 1971, Garde and Ranga Raju 1985), a 
satisfactory analysis from the point of view of basic mechanics is still out of reach. 
The present understanding of sediment transport is highly empirical in nature. In 
addition, numerical solution of three-dimensional equations of motion for water and 
sediment to simulate the actual movement of sediment is quit complicated, and is 
not necessary from the engineering point of view. Therefore, like most of the existing 
models, aggradation and degradation model presented in this study is based on the 
following considerations. 

1. We are only interested in computing the bed level changes and, not in com- 
puting the velocity of sediment particles. 

2. Sediment transport capacity can be satisfactorily quantified by a functional 
relation with average flow parameters. 
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With the above assumptions, a sediment continuity equation and an equation for the 
sediment transport capacity as a function of flow completely describe the sediment 
flow. Because sediment discharge is related to average flow parameters, it suffices 
to represent the water flow by cross sectional averaged quantities. 

As already mentioned in chapter - 1, it appears that quasi-steady flow uncoupled 
models do give satisfactory results in most of the cases of practical interest (Cui et 
al. 1996). The numerical model presented in this study belongs to that category. 
Quasi-steady flow models are based on the fact that, in general, the free surface 
wave celeraties are far greater than the celerity of bottom perturbation (Gunge et 
al. 1980). Therefore, it is possible to study bed perturbations while neglecting the 
time derivative terms in the water flow equations. This means that the water flow 
may be represented by a single steady state equation such as the energy equation. 

The other major assumption used in the present model is that the sediment 
is of uniform size. The size characteristic of the sediment is represented by the 
median diameter, d 50 . This assumption is not valid in cases where the sediment is 
well graded. During degradation, finer particles are transported in preference to 
larger particle. A stage may be reached when the bed has only sufficiently large 
particles which can not be moved by the flow. These large particles on the bed 
surface provide a protective layer for the underlying finer particles which would 
otherwise be transported. This phenomenon is known as armouring. Since the 
proposed model considers only a uniform sediment size, it can not be used where 
degradation is controlled by armouring. This aspect is being looked into a parallel 
study at IIT Kanpur. It should also be noted here that, like most of the currently 
available models, only bed load transport is considered and the effect of suspended 
sediment load is assumed to be negligible. 

In height of the above discussion, the essential feature of the aggradation and 
degradation models is the solution of four coupled relations : 
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1. Energy equation for water flow 

2. A friction-slope equation 

3. Equation for sediment discharge and , 

4. Equation of continuity for sediment. 

The following sections present details of these formulation aspects of the model. 

2.1 Equation for Water Flow 

Equation for water flow is based on the following assumptions : 

1. the pressure distribution in a cross-section is hydrostatic, 

2. the channel bottom slope is small, 

3. the energy correction factor, a = 1.0, 

4. all the channels in the system are prismatic to start with, and remain more or 
less prismatic even after aggradation or degradation and, 

5. the flow is steady. 

The steady state energy equation for a channel reach, as shown schematically in Fig. 
2.1, can be expressed as ( Chaudhry 1993, Subramanya 1997). 


dh % JZ S l 

dx i - 

1 gA? 


( 2 . 1 ) 


Where h — flow depth (m), Z = bed elevation with respect to a specified datum 
(m), Q = discharge (m 3 / sec), T = water surface width (m), A = cross-sectional 
area (m 2 ), g = acceleration due to gravity (m/s 2 ), Sf =the friction slope and, x = 
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distance along the channel (m). A general irregular cross- section as shown in Fig 
2.1b is considered in the present study. The channel width, B at any height y above 
the bed is assumed to be given as 


B = B 0 + C iy c > + C 3 y c < (2.2) 

where, Bo = the channel bottom width, Ci, C 2 , C 3 and C 4 are coefficients which 
are determined by the channel shape. The initial values of Ci, C 2 , C 3 and C 4 any 
for channel section at any distance may be obtained using the field survey data. 
Curve fitting methods may be used for this purpose. These values are continuously 
updated in the numerical model as the aggradation or the degradation of the section 
occurs. The top width, T is now given by 


T = B 0 + C x h n + C 3 h c 4 


The cross-sectional area A can be expressed as 


A = B 0 h + + 

1 + C-2 

The wetted perimeter P can be expressed as 


1 , 1 + 0 . 

1 + C, 


(2.3) 


(2.4) 


P = Bo + jf* [ACi C 2 )V 2C ’- 2 > + 1 + + l] dy (2.5) 

It should be noted here that C\ and C 3 are taken as zero in the case of rectangular 
channels. C 2 and C 4 are taken equal to one in the case of unsymmetrical trape- 
zoidal channels. The channel shape obtained using Eq. (2.2) is a cross between a 
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trapezoidal and an exponential channel. Eq. (2.2) and thus, the present numerical 
model cannot be applied for the case of channels with flood plains. An attempt is 
made in the present study to include the compound channel shapes also. However, 
the work remained incomplete due to lack of time. 

2.2 Sediment Continuity Equation 

The aggradation or the degradation in a channel reach is calculated by applying 
the sediment continuity equation between the two bounding sections. The sediment 
continuity equation is derived by applying the principle of conservation of mass to 
a differential control volume surrounding the channel bed. It should be noted here 
that only the bed load transport is considered in this study. The sediment continuity 
equation is expressed as (Raudkivi 1990) 

+ (2 - 6) 

where A s = the cross-sectional area surrounding the channel bed through which the 
bed load transport takes place, Q s = total bed load discharge through the cross- 
section, P = porosity and t = time. For a rectangular channel, Eq. (2.6) reduces to 


h n \ dZ + dq ° 

(1 “ p) aF + & 


where, q 3 = bed load transport per unit width. 


(2.7) 
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2.3 Sediment-Discharge and Friction-slope Pre- 
dictors 

There are no analytical equations which relate the friction-slope, 5/ and the sediment 
discharge Q s to flow and the sediment characteristics. Therefore, empirical equations 
are generally used. Several equations are available for this purpose. A good review 
of these equations can be found in the literature (Graf 1971, Vanoni 1975). Jansen et 
al. (1979) point out that many of the sediment discharge formulae can be expressed 
in a functional form as 


X = f(Y) (2.8) 

Where, A is a transport parameter which depends on the sediment discharge and 
the grain properties. Y is a flow parameter which depends on the flow properties 
and the representative grain diameter. The numerical model presented in this study 
is structured in such a way that any sediment discharge relationship can be easily 
considered. At the present stage of development, three alternative procedures are 
included in the numerical code, which are described as follows. 

t 

2.3.1 Power Law Model 

In this model, the sediment discharge is estimated by an empirical power function 
of the flow velocity. The unit sediment discharge, q 3 is given by 

«. = «(!)* (2-9) 

where a and b are empirical constants whose values depend upon the sediment 
properties. It should be noted here that Eq. (2.9) has been used earlier by many 
researchers (Soni et al. 1980, Zhang and Kahawita 1987, Bhallamudi and Chaudhry 
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1991). In the power law model, the friction-slope, 3/ is estimated using the Man- 
ning’s equation. 



Q 2 n 2 P a 

~a¥~ 


where, n = Manning’s roughness coefficient. 


( 2 . 10 ) 


2.3.2 Meyer— Peter— Miiller equation 

Meyer-Peter-Mfiller formula is very popular in Europe. It was developed based on 
experiments with coarse sand and gravel, of different relative densities. Most of the 
data upon which the formula is based were obtained in flow conditions with little 
or no suspended load. This equation can be expressed as (Chien 1954) 


$ 


r' 


(Ts Tta) ^50 


0.047 


1.5 


( 2 . 11 ) 


Qb j 'Yw ( 1 

Is) L(7s - Iw) d 5 0 \gd 5 Q 3 

where, r' = particle shear stress, d 5 o = median diameter of sediment particles, -y w = 
unit wt. of water, 7 S = unit wt. of sediment particles and qt, = sediment transport 
rate in weight per unit time per unit width. It should noted that Eq. (2.11) and 
Eq. (2.12) are in English units. As a first approximation, the particle shear stress, 
r' can be taken as bed shear stress rj which is given as 



0.5 


( 2 . 12 ) 


n = pg(^)s, (2.13) 

where, p = density of water. The friction -slope, 5/ in Eq. (2.13) may be estimated 
using Eq. (2.10) 
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2.3.3 Iowa Model 


Karim and Kennedy (1981) of the University of Iowa, developed a Total Load 
Transport Model(TLTM) for the sediment discharge and friction factor predictors. 
TLTM is derived by using non linear multiple regression analysis to a large body of 
laboratory and field data. This formulation takes in to account the well known fact 
that the friction-slope in alluvial channels is heavily dependent on the sediment 
discharge. The following sediment discharge and friction-slope predictors from 
TLTM are adapted as one of the alternatives in the present study. 

Sediment discharge predictor: 


log 


Qs 

.\M S ~ 


-2.2786 + 2.9719 log Vi 

+ 1.06 log Vi logU 6 + 0.2989 logU 2 logU 6 (2.14) 


Friction slope predictor: 


log 


Q 

A \Jg (s — 1) dso 


0.9045 + 0.1665 logU 7 

+ 0.0831 log V 4 logU 5 logVr + 0.2166 logV* logo's 
- 0.0411 log V 2 logU 3 log V 4 (2.15) 


where, 


Vi = 

Q 

(2.16) 

A\Jg(s — 1 ) d 50 

v 2 = 

h 

d 50 

(2.17) 

Vi = 

S f * 10 3 

(2.18) 

V, = 

jghS, 

U) 

(2.19) 
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( 2 . 20 ) 


V 5 




v 7 


to d^o 

V 

\J 9 h> Sf u+ c 
\Jg{ s - 1)^50 

<h 

y/g(s - l)4o 3 


( 2 . 21 ) 

( 2 . 22 ) 


In Eqs. (2.14)- (2.22), u> =fall velocity of sediment particles, v = kinematic viscosity 
of water, s = specific gravity of sediment particles and u* c = critical shear velocity 
obtained from the Sheilds diagram. Following explicit equation for u* c (Raudkivi 
1990) is used in the present study 


d* = 2.15 R e » 

Re. 

< 1 

(2.23) 

d, = 2.5 (R e ,) a8 

1 < R e » 

< 10 

(2.24) 

d. = 3.8 (Re*) 0 ' 625 

Re* 

> 10 

(2.25) 

In Eqs. (2.23-2.25), d* = dimensionless particle 

size = [ (7s - 

- 7)Wso 3 A' 2 ]3, and 


Re* = particle Reynolds number = yJg^Sj * ^ L 


In all the three prior discussed models, equations are available only for deter- 
mining the unit sediment discharge. Therefore, these models can be readily applied 
in case of rectangular channels and wide channels. It should noted here that no 
predictive equations are available for estimating total sediment dischargers) in 
narrow irregular channels. Q s in these channels is assumed to be given by the 
following equation. 


Qs = T q s C 5 


(2.26) 


where, C5 is an empirical coefficient which depends upon the shape of the cross- 
section. C 5 is taken equal to one in all the simulations reported herein. 
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Chapter 3 


Computational Procedure 


Several numerical techniques are adopted to solve the governing equations of quasi- 
steady flow, uncoupled model for application to tree-type alluvial channel networks. 
The details of these procedures are discussed in this chapter. The overall algorithmic 
structure of the model is also presented. 


3.1 Determination of Bq, C\, C2, C3 and C4 

The coefficients Bo, Ci,C 2 , C 3 and C 4 in Eq. (2.3) are determined from the known 
values of widths B(y) at some discrete points along the x-axis. Linear regression 
analysis is adopted for this purpose. Width, Bi(y) is divided into two parts B r and 
Bi where, B r — the distance between the right bank and the y- axis while Bi = the 
distance between the left bank and the y- axis. The curve for the left bank can now 
be expressed in a functional form as 

Bi{y)-B l (0) = C 1 y°* (3.1) 

Equation (3.1) may be transformed to a linear relationship as follows 

x , = C , 1 + C 2 y' ^ (3.2) 
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where, x' = \n[Bi(y) — Bi( 0)], y' — ln(y), C[ = ln(Ci). C{ and C 2 can be 
determined from the known values of x r at some discrete points along the x- axis 
using the least squares linear regression (Press et al. 1986). The equations are 


and 


<?2 


Mjf s' - gfei yjarj 

Mff)* - vf 


(3.3) 


C; = x' - C 2 y’ (3.4) 

where M — the total number of points at which £;(y) is known, and the overbar 
indicates the mean values. Value of C\ can be easily determined from C[. The same 
procedure is adopted for determining C 3 and C 4 from the known values of B r (y). 
It should be noted that the values of B 0 , C\, C 2 , C 3 and C 4 are updated at the 
end of each computational time step because the aggradation or degradation at the 
cross-section may alter the values of Bi(y ) and B r (y). 


3.2 Computation of Wetted Perimeter, P 

The cross -sectional area, A corresponding to a known value of flow depth, h can be 
easily computed using Eq.(2.4). However the computation of the wetted perimeter, 
P using Eq. (2.5) is not trivial. Five point Gaussian quadrature integration as 
described below is used for this purpose. Equation (2.5) may be recast as 

P = B 0 + [ f(y) dy (3.5) 

J 0 

where 


f(y) = + + v'l + (C' 3 C 4 ) 2 y 2C <- 2 (3.6) 


24 



k 

1 

2 

3 

4 

5 


0.236926885 

0.478628670 

0.568888889 

0.478628670 

0.236926885 

Y k 

-0.906179846 

-0.538469310 

0.0 

0.538469310 

0.906179846 


Table 3.1: W k and Y k for the Gaussian quadrature 


Eq. (3.5) can noe be approximated by the five point Gaussian quadrature as (Press 
et al. 1986) 

■ P = \ £ W k f(.y' k ) (3.7) 

Z k = 1 

where, W* = value of the weight corresponding to the k’th point, and y' k = normal- 
ized value of y at the k’th point. 

y'i = |(1 + n) (3.8) 

W k and Y k values for k = 1 to 5 are given in the table 3.1 


3.3 Computation of Sediment— discharge 
and Friction-slope 

Computation of sediment-discharge, q s and friction-slope , 5/ for the known flow 
and sediment characteristics is trivial in the case of Power law and the Meyer-Peter- 
muller models. However, this is not so in the case of TLTM proposed by Karim 
and Kennedy (1982). Equations (2.14) and (2.15) are solved using an inefficient 
first-order iteration technique in the IALLUVUIAL model devopled by Karim and 
Kennedy (1982). In this study, the Newton-Raphson technique is used to solve for 
the q s and Sf in the TLTM model. 
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A single equation for the unknown Sf is first derived by eliminating g 5 from Eqs. 

(2.14) and (2.15). Denoting 


<Jg(s-l)d 5Q = ai 

(3.9) 

'sfgh = a 2 

(3.10) 

loe{ Aa? = ° 3 

(3.11) 

log '(-T-) = 04 

«50 

log( , ) = x i 

a l &50 

a 2 

(3.12) 

(3.13) 

(3.14) 


Eq. (2.14) can be written as 


xi = [-2.2786 + 2.9719 a 3 - (1.06 a 3 + 0.2989 a 4 ) log(ai)] 

+ (1.06a 3 + 0.2989 a 4 ) log(x 2 - u tc ) (3.15) 


Denoting 


-2.2786 + 2.9719 a 3 - (1.06 o 3 + 0.2989 a^log^) = a 5 (3.16) 


and 


1.06 a 3 + 0.2989o 4 = a 6 


(3.17) 


Eq. (3.15) can be written as 


Xi = a 5 + a 6 log (a ; 2 - u* c ) 


(3.18) 


Denoting 


log 


(^) = 


a 7 


(3.19) 
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, ( 1000 \ , x 
° S \aT*aJ = 08 ( 3 - 2 °) 

0.1665 — 0.0831 a 7 log(w) = ag (3.21) 

0.2166 a 7 - 0.0411o 8 + 0.0822 log(co) = o 10 (3.22) 

0.0831 07 = On (3.23) 

- 0.0822 o 4 = 012 (3.24) 

0.9045 - 0.2166 a 7 log(co) + 0.0411 a 4 a 8 log(a;) - a 3 = o 13 . (3.25) 

and substituting Eqs. (3.9)-(3.14) and (3.19)-(3.25) in Eq. (2.15) one can obtain 

agoq + a 10 log(x 2 ) + on xi log(x 2 ) + a 12 {log(x 2 )} 2 + a i3 = 0 (3.26) 


Denoting 

Og Os + 0 43 = bi 

09 a$ = b 2 

O10 + 05 an = 63 

Oil a 6 — bi 

and substituting Eq. (3.18) in Eq. (3.26) one can obtain 

/(x 2 ) = 61 + 6 2 log(x 2 - u, c ) + 63 log(x 2 ) 

+ bi log(x 2 ) log(x 2 - u +c ) + ai 2 (log(x 2 )) 2 = 0 (3.31) 

Coefficients 61, i> 2 , b 3 , bi and oi 2 in Eq. (3.31) are functions of known flow and 
sediment characteristics, h, Q, d 50 , v and to. The critical shear stress, u* c is also 
determined explicitly using the median diameter of the sediment particles, d 50 . The 
only unknown in Eq. (3.31) is x 2 , and it can be solved for using the Newton-Raphson 
technique. In this method, a value is first assumed for the x 2 , and is then improved 
iteratively using the following equation, 


(3.27) 

(3.28) 

(3.29) 

(3.30) 



where, 


X 2 


x 2 old - 


f(x 2) 

I'M 


(3.32) 


f'M 


df(x 2 ) 

dx 2 

b 2 , 63 , &4log(x 2 - M*c) 

-j- -J- 

3?2 ^*c ^2 #2 


^4 X2 
X2 U* c 


+ 2 aw 


Iog(x 2 ) 

2:2 


(3.33) 


A very small value for the error tolerance is used while solving the Eq. (3.31)by the 
Newton-Raphson technique. The friction-slope, Sf is computed using Eq. (3.14) 
once the value of x 2 is determined. Subsequently, q s can be determined using Eqs. 
(3.18) and (3.13). It should be noted here that this minor modification to the way 
the TLTK1 is applied results in significant savings in the total computational time 
because Sf needs to be computed at least four times at each section during each 
computational time period. 


3.4 Solution Of Gradually Varied Flow Equation 

Equation (2.1) representing the steady gradually varied flow in open-channels is 
used to determine the flow depth at any section ‘j’ in channel ‘i’ once the flow 
conditions at the downstream section ‘j+1’ in the same channel are known. Fourth- 
order Runge-Kutta method (Chaudhry 1993, Subramanya 1997) is used to solve Eq. 
(2.1) numerically. 

At 

Kj = h w - (SL n + 2 SL 2 + 2 SL Z + 5L 4 ) (3.34) 
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where, the slopes SL X , SL2, SL$ and 5L 4 are determined using the following pro- 
cedure. 

1. For the known discharge, Qi and the flow depth at section i+1, compute 
the flow area, A, top width, T, and the wetted perimeter, P using the geometric 
properties of the cross-section (ij+1). Eq. (2.3), (2.4) and the procedure given 
in section 3.2 are used for this purpose, respectively. 

2. Sf at section (i,j + 1) is computed using either the Manning’s equation or the 
procedure described in section 3.3 for the TLTM model. 

3. Above values of A, T and Sf are used in Eq. (2.1) to compute the SL X = — 
value. 

4. A predicted flow depth, bl is determined as 

A nr 

h! = hi J+1 - — (SL X ) (3.35) 

where, Aa; = the distance between the section (i,j) and (ij+1). 

5. Using the h! value as determined from Eq. (3.35) and the geometric properties 
of the section midway between (i,j) and (i,j+l), steps 1-3 are repeated to 
determine the value of SL^. 

6. A new predicted flow depth, h" is determined as 

Ar 

h" = h iJ+l - — (SLi) (3.36) 

7. Steps (l)-(3) are repeated using h" and the section properties at (i, J + |) to 
determine the value of SL 3 . 

8. A predicted flow depth, h"' is determined as 
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(3.37) 


h!" — hij+i — Ax (SL3) 

9. steps 1-3 are repeated using h'" and the section properties at (i j) to determine 
the value of SX4 

Runge-Kutta method is applied in a marching fashion from the node(i,n,) to the 
node (i,l) on the channel i. To start the computations, the discharge in the channel, 
Qi and the flow depth at the last node of the channel, h itTli should be specified. h^ ni 
is known either from the specified downstream boundary condition for a pendant 
channel or from the previous computations. 

3.5 Computation Of GVF Profiles in Channel Net- 
works 

A large part of the computational expenditure in the application of the aggradation 
or degradation model to alluvial channels is incurred by the water surface profile 
computations. This is especially true in the case of large channel networks. There- 
fore, the major emphasis in the present study is to perform these computations 
as efficiently as possible. Application of the fourth-order Runge-kutta method as 
described in section 3.4 is straight forward in the case of single channels. However, 
this is not so in the case of channel networks. In this study two separate computer 
codes are developed for (1) the converging tree-type channel networks and, (2) 
the diverging tree-type channel networks, respectively, by exploiting the special 
structures inherent in these problems. 
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Node number 



3.5.1 Converging Tree- type Channel Networks 

In a typical converging tree- type channel network (Fig. 3.1), flows from several 
tributaries join together to form the flow in the main channel. This is typical of 
river systems in non-deltaic regions. For the sake of simplicity, it is assumed that 
(1) discharge from the system leaves through one single channel at the downstream 
end, (2) there are no hydraulic structures within the system, (3) only two channels 
feed into single channel at a junction and, (4) the flow is subcritical throughout the 
system. 

Assumptions (1) and (3) apply to most of the natural channel systems while the 
assumption (2) can be easily relaxed to include intermediate control structures. The 
occurrance of the supercritical flow in channels is rare and requires special treatment. 
Analysis of such a system is beyond the scope of the present study. The problem of 
GVF computation in converging tree-type channel systems can be defined as 
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“given the discharges in all the root channels( Qi, Q 2 , Qa and Q 5 in Fig. 3.1) 
and the flow depth at the pendant node(/ig in Fig. 3.1), compute the water surface 
profile in the system” . 

The above problem can be solved by the successive application of the procedure 
for single channel, termed as the Initial Value Problem(IVP) hereafter, and by 
applying the compatibility conditions at a junction. The compatibility conditions 
at a junction are given by 


Qi + Qj — Qk 


(3.38) 


and 


hi,ni “I - Z i,m — hj'Tij ~t~ Zj t nj — hk t i + (3.39) 

where, the subscripts i and j refer to the incoming channels at a junction and the 
the subscript k refers to the outgoing channel, n* and rij refer to the last nodes on 
channels i and j, respectively. Equation (3.38) refers to the mass balance equation 
and Eq. (3.39) refers to the energy balance equation. It should be noted here that 
the energy losses at a junction and the velocity heads are neglected in Eq. (3.39). 
However, these assumptions can be easily relaxed if required. Solution for the IVP 
is applied successively for the converging networks as described below. Channel 
network shown in Fig. 3.1 is considered for the purpose of illustration. 

1. Discharge in channel [7] is equal to the total discharge in all the root channels 
i.e.,<3i + Q 2 + Qa + Qs- Flow depth at the node 8 is known either from the 
specifed boundary condition or from the rating curve at that section. 

2. IVP is solved for channel [7] using the method described in section 3.4 and, 
Ha/ 7, the flow depth at node 4 in channel [7] is determined. 
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3. Equation (3.39) is used to determined the flow depths h 4 / 3 and h 4 / 6 in channels 
[3] and [6] respectively. In channel [6] i.e., Q 6 is equal to the total discharge 
in all the root channels feeding to it i.e.,Q 6 = Q 4 + Q 5 . Discharge in channel 
[3], Q3 — Q7 — Qs~ 

4. IVP is now applied to channels [3] and [6] to determine h 2 / 3 and h 6 / 6 , respec- 
tively. 

5. Step-3 is applied to nodes 2 and 6 to determine (h 2 / 2 , Q 2 ); (hi/ 2 , Qi); (h 6 / 4 , Q 4 ); 
and (h 6 / 5 , Q 5 ) 

6. IVP can now be applied to channels [1],[2],[4] and [5] to complete the water 
surface profile computation in the whole network. 

The above discussed procedure is used to develop a computer code for determining 
water surface profiles in general, converging tree-type channel networks, given (1) the 
inflow discharges in all the root channels, (2) the flow depth at the pendant node, (3) 
bed elevation variation in all channels and, (4) variation in channel section properties 
in all channels. This is used as a subroutine in the aggradation or degradation model 
for the converging channel networks. The subroutine presented in this section itself 
calls the procedures presented in sections 3.2, 3.3 and 3.4 as subroutines. 

3.5.2 Diverging Tree-type Channels Networks 

In a typical diverging tree-type channel network (Fig. 3.2), flow from a root channel 
gets distributed to several branches. Such networks are common in (1) irrigation 
systems and, (2) river systems in deltaic region. For the case of simplicity ,it is 
assumed that (1) discharge from a single source enters the channel at upstream 
end, (2) there are no hydraulic structures within the system, (3) there are only two 
branches to a parent channel and, (4) the flow is subcritical through out the system. 
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However, the principles presented in this study can be easily extended to networks 
with multiple source and intermediate control structures. They can also be extended 
to cases where there are more than two branches to a parent channel. The problem 
of GVF computation in diverging tree-type channel systems can be defined as 

“Given the discharge in the root channel ( Q 2 in Fig. 3.2) and the controlling 
depths at all the pendant nodes (h 4) h 5 , h 8 , h 9 , tin and h 12 in Fig. 3.2), compute the 
water surface profile in the system. Variation in the bed level and the channel shape 
parameters are also specified for all the channels.” Successive application of the 
IVP as described in section 3.5.1 for this case would involve a trial and procedure 
because the discharge in the pendant channels is not known apriori. In this study, 
an efficient algorithm developed by Naidu et al. (1997) is used to compute the 
GVF profile in a diverging tree-type channel network. Details of this algorithm are 
presented elsewhere (Naidu 1995 Naidu et al. 1997) and only a brief description 
is presented here for the sake of completeness. The network shown in Fig. 3.2 is 
considered for the purpose of explanation. 

The first step in the preposed method is to assume the discharge at a pendant 
node, say node 12. With this, the discharge and the controlling depth at the 
downstream end for channel [12] are known. Ttherefore, Eq. (1) can be solved for 
channel [12] as an “Initial Value Problem” (IVP) by marching the computations from 
node 12 to node 10. The Fourth-Order Runge-Kutta method is used for this purpose. 
At junction node 10, Eq. (3.39) can be applied to determine the upstream depth 
in channel [11] since the upstream depth in channel [12] is known from previous 
computations. With this, both upstream and down stream depths, i.e., at nodes 
10 and 11 respectively are known for channel [11]. However, the discharge is not 
known. Therefore, Eq. (1) is solved for channel [11] as a “Boundary Value Problem” 
(BVP) for a single channel using the ’’shooting method” (Press et al.1986). In this 
method, the BVP is solved as an IVP with iterations.This computation determines 
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the discharge in channel [11] as a part of its solution. The continuty equation, 
Eq. (3.38) and the energy equation, Eq. (3.39) can now be applied at node 10 for 
determining the discharge and downstream depth for channel[10]. Computations can 
then be marched in channel [10] from node 10 to node 6 as an IVP, Eq. (3.39) can 
be applied to determine the u/s depth of channel [7]. However, neither the discharge 
nor the downstream depth in channel [7] are known. Therefore, computations at 
junction node 6 are different from the computations at the junction node 10. In t his 
case, the group of channels [7], [8] and [9] are considered simultaneously. For the this 
subsystem, flow depths at nodes 6,8 and 9 are known. Therefore, it can be solved as 
a BVP for a group of channels. The computations are carried out as follows. 

• Assume discharge in channel [8] 

• Apply IVP in channel [8] 

• Apply junction conditions at node 7 

• Apply BVP in channel [9] 

• Apply IVP in channel [7] 

• If the computed depth at node 6 is same as the previously determined flow 
depth at this node, the assumed discharge in channel [8] is correct. Otherwise, 
it is updated using the ’’Shooting Method”. 

The above computations give the discharge in channel [7]. Equations (3.38) and 
(3.39) can now be applied at node 6 to determine the discharge and downstream 
depth in channel [6]. An IVP can be solved now to march the computations in 
channel [6] from node 6 to node 2. It can be seen from Fig. 3.2 that a BVP for a 
group of channels has to be solved at node 2 in order to determine the discharge 
in channel [3]. If this computed discharge in channel [2] is sum of the discharges 
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in channels [6] and [3]. If this computed discharge in channel [2] is same as the 
prescribed discharge at the root, the sssumed value of discharge in channel [2] is 
corerect. Otherwise, it is Updated using, once again the the shooting method and the 
computations are repeated till convergence. The flow chart for the above procedure 
with reference to the example network (Fig. 3.2) is shown in Fig. 3.3 



Figure 3.2: Typical diverging tree-type network 
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Figure 3.3: Flow chart for GVF Profile in diverging tree-type network 
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It is obvious from the above procedure that different “Paths of Marching” for the 
solution may be obtained depending on the pendant node at which the computations 
are started. The path of marching also depends on node chosen for starting the iter- 
ations of “Boundary Value Problem for Group of Channels, BVPGC” . The proposed 
algorithm results in high efficiency when calculations march along a particular path 
in the network. For example, consider the case of starting the computations at the 
pendant node 5 for the network shown in the Fig. 3.2. This choice would result in 
a six fold increase in the computational time as compared to the previous choice of 
starting the computations at node 12. Naidu et.al. (1997) presented an algorithm 
for the best path of marching using some simple graph theoritic principles. 

The computer code developed by Naidu et al. (1997) is used in this study to 
compute water surface profiles in diverging tree type channel networks. The original 
code developed by Naidu et al. (1997) used a uniform slope for each channel and 
considered only cross-sections of trapezoidal shape. This code is modified by the 
author to include variable bed slope and more general cross-sectional shapes as 
given by Eq. (2.2). The modified code is used as a subroutine in the aggradation/ 
degradation model for diverging tree-type channel networks. 


3.6 Computation of Aggradation/Degradation 

The aggradation/degradation at a channel section is computed by numerically solv- 
ing the sediment continuity Eq. (2.6). While solving the Eq. (2.6), prior computed 
sediment discharge values from the water flow module are used. Backward spatial 
and backward time finite- differences are used to discretize Eq. (2.6). The change 
in the sediment flow area at any any section *j’ in channel V, A A SiJ . is given by 

“ l - p fa; ~ G«ij] ( 3 - 40 ) 
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or, for a wide channel, 

& Z ij = l _p Arc ~ 9si ->] ( 3>41 ) 

where, At = computational time step. Equations (3.40) and (3.41) should not be ap- 
plied at the upstream nodes of the root channels and at the junction nodes. Bound- 
ary conditions are applied at the upstream nodes of the root channels. Junction 
nodes are treated in a special manner to ensure sediment mass balance condition. 

3.6.1 Upstream Boundary Condition 

Inflow sediment discharge hydrographs are specified as the boundary conditions at 
all the upstream nodes i.e. j=l in all the root channels. i is determined using 
the following equation. 

A ^,< = K. - ««,] ( 3 - 42 ) 

where Q Si 0 refers to the sediment inflow at time t+At for channel ‘i’ and is obtained 
from specified inflow sediment discharge hydrograph. A similar procedure was used 
by Bhallamudi and Chaudhry (1992) and satisfactory results were obtained for the 
case of aggradation due to sediment overloading. 

3.6.2 Aggradation/Degradation of Junction nodes 

The procedure for determining the aggradation/degradation at a junction node is 
slightly different from Eq. (3.40) for an interior node. Refering to Fig. 3.4, 


AZjun 


Qs tl,n f i — 1 Qs>3,2 

B'n Axil + B ' i2 A x i2 + B' i 3 Ax i3 


(3.43) 
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Figure 3.4: Schematic layout of a junction 


where, A Zj un = change in the bed level at the junction in time At. Subscripts 
(il,riii — 1) and (*2, n i2 — 1) refer to penultimate nodes in channel il and i2, 
respectively. Subscript (i3,2) refers tc the second node in channel i3. Subscripts il,i2 
and i3 refer to the channels il,i2 and i3 respectively. B' refers to an “effective bottom 
width” for the channel. This concept is necessary to avoid singularity problem in 
case of channels of zero bottom width such as, the triangular channels. At the 
present stage of code development, B' is taken as the bottom width B 0 . Numerical 
simulations are carried out only for those cases where all the channels have finite 
bottom width. It can be seen that Eq. (3.43) is an expression for sediment mass 
balance condition for a control volume sarrounding the nodes (il, nn — 1), (i2, n i2 — 1) 
and (i3, 2). An inherent assumption in Eq. (3.43) is that the bed aggrades/ degrades 
by equal amounts on all the three sides of a channel junction. This assumption is 
questionable in case of junctions with free-overfalls and needs further study. 

It should also be noted here that the aggradation/degradation at the last node of 
any pendant channel is computed also using Eq. (3.40). The inherent assumption in 
this is that the outflow sediment discharge from a pendant node is equal to sediment- 
discharge capacity at the last node of the channel. This assumption is generally valid 
for the test cases presented in the study. However, it may not be strictly valid in 
some other cases. For example, if the downstream control structure is a spillway, it 
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may allow water to flow over, but reduces sediment outflow to zero. 

3.6.3 Change in Cross-Section 

Aggradation/degradation at a cross-section will result in change in the bed level as 
well as the cross-sectional shape. For wide channels, these changes are completely 
described by Eq. (3.41) since there would not be any appricable change in the 
bank configuration. For non-rectangular channels, however, the value of AA Sij . at 
any cross-section should be used to update the bed levels and the channel shape 
parameters. 

Aggradation: In the case of aggradation, it is assumed that the banks would not 
be able to support the deposited material and deposition during the time At occurs 
completely on the bed (Figure 3.5a). Change in the bed level, A Z is computed 
inversely from the following equation. 

AA S = Bo AZ + -^-(AZ) l+C2 -^-(AZ) 1+c * (3.44) 

1 TO2 1 ~r O4 

Updated value of Bq is determined using the equation 

B™ w = B° 0 ld + Ci(AZ) C2 + C z (AZ) Ci (3.45) 

Values of Bi and B r are accordingly changed and new values of Ci, C 2 , C 3 and C 4 
are determined using Eqs. (3.3) and (3.4). Slight deviations in the bank profile may 
be observed because of the statistical procedure adopted in determining Ci, C 2 , C 3 
and C4. 

Degradation: Degradation is assumed to occur due to erosion from the entire 
wetted perimeter. Two alternative procedures are adopted to compute the new cross- 
section due to degradation. In the first alternative (Fig. 3.5b), a parallel degradation 
is assumed, i.e., the entire cross-section expands uniformly by an amount AZ. AZ 
in this case is given by 
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(a) Aggradation 



(b) Degradation : Alternative - 1 



\ / 
v / 


(C) Degradation : Alternative - 2 

Figure 3.5: Aggradation/Degradation in a non-rectangular channel 
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AZ = 


A A, 

P 


(3.46) 


Bq, Bi(y ) and B r (y) values are oocordingly changed and new values of C x , C 2 , C 3 
and Ci are determined using Eqs. (3.3) and (3.4). 

In the second alternative, the degradation at any point on the banks or the bed 
depends on the local shear stress. The local shear stress r[ is assumed to be given 
by (Chang and Hill 1982) 


r[ = pgh! Sfcosd (3.47) 

where, h! = local flow depth and 9 = local transverse slope of the bank. 

The decrease in area will be distributed such that the local scour depth ( Z s ) is 
proportional to the local shear stress. The local scour depth can be obtained by 
(Chang and Hill 1976) 


^ = 


- r CT 

£r fa' ~ T cr ] Ay 


AA 


(3.48) 


where, r cr — critical shear stress, Z s = local scour depth Ay = the horizontal 
distance between 2 points and T = top width of the channel. If the shear stress is 
negitive then it is set to zero. The new values of C\, C 2 , C 3 and C 4 are determined 
using Eqs. (3.3) and (3.4). 


3.7 Quasi-Steady Flow Uncoupled Algorithm 

The principles presented in sections 3.1 through 3.6 are used to develop quasi-steady 
flow uncoupled models for aggradation/degradation in alluvial channel networks, the 
flow chart of the algorithm is presented in Fig. 3.6. 
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Subroutine for 

(1) Power low or 

(2) MPM model or 

(3) TLTM 


Figure 3.6: Flow chart for the quasi-steady flow model 
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3.8 Closure 


In this chapter, details of quasi-steady flow uncoupled models for computing aggra- 
dation and degradation in tree-type alluvial channel networks are presented. It 
is envisaged that these models will be eventually used for simulating bed level 
changes in networks as large as the Ganga river system, and for periods as long 
as 20 years. Therefore, every effort is made is to reduce the computational cost per 
one computational time period while formulating these models. 
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Chapter 4 


Results And Discussion 


In this chapter, the quasi-steady flow, uncoupled model is verified by comparing 
the present numerical results with (1) observed data (Soni et al. 1980) and ear- 
lier numerical results (Bhallamudi and Chaudhry 1991) obtained using a complete 
unsteady, coupled model. The applicability of the numerical model to converging 
and diverging tree-type alluvial channel networks is demonstrated by simulating 
(1) aggradation due to sediment overloading, and (2) degradation due to sediment 
depletion. Results are presented for both rectangular channels and nonrectangular 
channels. 

4.1 Verification of The Model 

The present numerical model was used to simulate the aggradation process observed 
by Soni et al. (1980) in a laboratory flume due to sediment over loading at the inflow 
section. Soni et al. (1980) used a rectangular flume with alluvial bed carrying a 
constant unit discharge, (fa, at a uniform flow depth of h 0 . The equilibrium between 
the water flow and the sediment discharge was distrubed by increasing the sediment 
discharge at the upstream end from its equlibrium value of q s0 to q 3 o + A q a . This 
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resulted in the aggradation of the flume bed so that the increased bed slope could 
transport the additional imposed sediment load that deposited on the bed. 

The experiments were conducted in a flume 0.2m wide and 30m long. The sand 
forming the bed and the injected meterial had a mean diameter, d 50 of 0.32mm. The 
values of empirical constants a and b in the sediment transport formula Eq. (2.9) 
as found from the uniform flow experiments were 1.45 x 10~ 3 and 5.0, respectively. 
Manning’s roughness coefficient ‘n’ was approximatly equal to 0.022 and the porosity 
‘p’, of the sediment bed layer was equal to 0.4. 

In the numerical simulation, the values of bed slope so = 0.0026, unit discharge, 
<?o = 0.032 m 2 /s, flow depth, ho = 0.075 m and, sediment load increment, Aq s — 
2.7q s0 . These data were taken from the run identified as test no. 2 by Bhallamudi 
and Chaudhry (1991). The initial bed elevations at all finite difference nodes as 
calculated from s 0 were specified as the initial conditions. Soni et al. (1980) 
conducted the experiments in a single channel. However, a three channel system 
as shown in Fig. 4.1 was considered for the application of the numerical model. 
The width of the pendant channel [3] was adjusted(B 3 =0.4m) such that the flow 
conditions in the two root channels [1] and [2] were almost same as the flow conditions 
obtained by Soni et al. (1980). Lengths of the channels were 30m each, with the 
total number of reaches in a channels being 30 (Ax = 1.0m). The computational 
time step, At was taken as 1.0 sec. The uniform flow depth was specified as the 
downstream boundary condition and sediment loads q 3 o + A q a were specified as the 
two upstrem boundary conditions for the root channels [1] and [2]. 

The numerical run was made for a total time, t=3000sec. The bed disturbance 
i.e., the aggradation profile, from the root channels did not reach the confluence 
point during this period of computation. Therefore, the aggradation characteristics 
in the two root channels in the numerical simulation were close to those in Soni et 
al. (1980). 
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A Z(x,t)/AZ(x,0) 



Figure 4.2: Variation of relative aggredation with 77, sediment overloading 
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Figure 4.2 shows the variation of the relative aggradation, i.e., A Z(x, t)/AZ(0, t ) 
with the- non-dimensional parameter, 77 = x/2y/Ktf, where K Q is the diffusion 
coefficient which is equal to bq s0 /3so{l ~ p)- Numerical results obtained using the 
present model are compared in Fig. 4.2 with the measured data (Soni et al. 1980) 
and with the numerical results obtained by BhaUamudi and Chaudhry (1991) using 
a complete unsteady coupled model. It can be observed from Fig. 4.2 that the non 
dimensional bed profiles obtained by the present quasi-steady flow model compared 
very well with the results obtained by Bhallamudi and Chaudhry (1991) who used 
a complete one dimensional partial differential equations for un-steady flow in a 
wide rectangular channel with deformable bed. This satisfactory matching indicates 
that the quasi-steady flow assumption and the uncoupled solution technique do not 
introduce any significant errors when simulating bed level changes. This aspect of 
modeling was confirmed also by Cui et al. (1996). However, the quasi-steady flow 
model may break down if the Froude number for the flow is close to unity. Figure 
4.2 is drawn for aggradation in root channel [ 1 ]. Identical results were obtained for 
the other root channel [ 2 ] as expected. The bed wave had not reached the junction 
during the period of computation. Therefore, the bed profile in the pendant channel 
[3] did not change as expected. 

4.2 Aggradation Due to Sediment overloading 

The objective of this work is to develop a numerical model for aggradation in genaral 
diverging and converging tree-type networks which considers genaral nonrectangular 
cross-sectional shapes and a sophisticated predictor for computing sediment trans- 
port capacity. To demonstrate its applicability, several test runs were made, results 
of which are discussed in this section and the following section . This section presents 
the results for aggradation due to sediment over loading. 
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4.2.1 Network With More Than One Junction 

In the previous section, simulations were made for a three channel system with 
one junction point. However, in a general network there can be more than one 
junction. Simulations for Scenarios (1) and (2) in this study were made for a five 
channel system of rectangular channels as shown in Fig. 4.3. The Power law model 
and the Iowa model for sediment discharge were used in Scenarios (1) and (2), 
respectively. The lengths of all the five channels were taken as 30m. An initial bed 
slope so=0.0026 was taken for all the five channels. The values of bed elevations as 
calculated from sq were assigned to each finite difference node as initial conditions. 
The constants a and b in Eq. (2.9) were taken as 1.45 x 10 -3 and 5.0 respectivly. 
The discharge per unit width in all the channels was taken equal to 0.032m 2 /s. 
Mannings’s roughness coefficient in all channels was 0.22. The width of each root 
channel was equal to 0.2m. Widths of other channels in the system were taken such 
that the unit discharge was constant through out the system. Sediment overloading 
was taken as 3.7 times the initial sediment carrying capacity at the upstream most 
node in all the root channels. Porosity of the bed layer was taken as 0.4. 
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The model was run for 3000sec with At=1.0sec and Az=1.0m. Figures 4.4 and 
4.5 show the initial bed profile and the final bed profile at t=3000 sec. for channels 
[1] and [2], respectively. 

The bed wave from the root channels had not reached any of the junction points. 
Therefore, there should not be any change in the bed profile in the channels other 
than the root channels i.e., channels [2] and [3] in Fig. 4.3. The bed profiles for the 
channel [2] shown in Fig. 4.5 confirms that the model gives consistent results. The 
bed profiles for all the root channels were identical to each other as expected. The 
bed profiles for channels [1] are presented in Fig. 4.4. It is obvious from this figure 
that the bed wave had not reached the junction point. 

In the simulation for the Scenario 2, the initial and other conditions were identical 
to those in Scenario 1, except the model for the sediment carrying capacity. The 
Iowa model was used to compute Sf and q s in this case. The mean size of the 
sediment particles, d5o=0.32mm. Figures 4.6 and 4.7 shows the initial bed profile 
and the final bed profile at t=3000sec in channel [1] and [2] respectivly. AS expected 
there was no aggradation in channel [2] Fig. (4.7) because the the bed wave had not 
reached the junction points. One can not compare the results obtained for Scenario’s 
1 and 2 because two different empirical models were used for computing the sediment 
carrying capacity. However, the constiency in the results can be checked. It can be 
observed from Figs. 4.4 and 4.6 that the aggraded more in the case of Scenario 
1. This is because the sediment carrying capacity and consequently the sediment 
overloading is more if the power law is used. Sediment carrying capacity at the 
upstream most section, q s0 is equal to 1.75 x 10 ~ 5 m 2 /s if Eq. 2.9 was used, while 
it is equal to 7.25 x 10~ 6 m 2 /s if Eqs. 2.14 and 2.15 were used. 
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Bed elevation 


t = 3000 sec. 
t=0 sec. 


Distance alogle channel^ 12 


Figure 4.6r Bed profile in Channel 1, Scenario 


t=3000 sec. 
t=0 sec. 


Distance along die channe?(m) 12 

Figure 4.7; Bed profile in Channel 2, Scenario. 




Figure 4.8: Bed profile in Channel 1, Scenario-3 


4.2.2 Aggradation Beyond The Junction points 

In Scenarios 1 and 2, it was observed that the bed wave had not reached the junction 
points. The simulation for Scenario 3 was performed to demonstrate the applicability 
of the model when the bed wave propagates beyond the junction point. The test 
conditions for the Scenario 3 were identical to those in Scenario 1 except that the 
lengths of the channels were 17.0m in this case. Also, the computations for the 
Scenario 3 were performed for a total time of 7000 sec. 

Figures 4.8 to 4.12 shows the initial bed profile and the bed profile after 7000 
seconds in channels [1],[2],[3] ,[4] and [5] respectly. It can be observed from Figs. 
4.8 to 4.12 that the aggradation had occured in all the channels. This indicates 
that the bed wave had moved beyond the junction points. No numerical difficulties 
were encountered while performing the computations. It can be observed from 
Fig. 4.9 that there is a sharp rise in the bed level change at the downstream end 
of the channel [2]. This is because the bed wave from the channel [5] arrived at 
the junction 2 earlier than the disturbance from the upstream side. The resulting 
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Figure 4.13: Bed profile in Channel 1, Scenario-4 


aggradation moved in an upstream direction in channel 2 toward junction 1. This 
occured simultationeously with the movement of aggradation in the downstream 
direction from junction 1 to 2. 

4.2.3 Effect Of Computational Time Step 

In the simulation for the Scenarios 1, 2 and 3 the computational time step, At 
was equal to 1.0 sec. Simulation for the Scenario 4 was performed to demonstrate 
that the results can be obtained with a larger At without introducing significant 
error. In the Scinario 4 the test conditions were identical to Scenario 3, except that 
At=60sec. Figures 4.13-4.17 show the initial bed profile, the bed profile after 7000 
sec with Af=1.0sec and, the bed profile after 7000sec with Af=60sec. It can be 
observed from these figures that even when At was taken as large as 60 times the 
previous value, the difference in the aggradation profiles was negligible. This fact 
is important from the point of view of applicability of this model to larger networks. 
The courant number for stability in the application of a quasi-steady flow model 
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.. ' Table 4.1: Cross-sectional shape parameters for Scenario-5 


channel no. 

Cross-sectional Shape Parameters 

B o(m) 

Ox 

<? 2 

C 3 

c 4 

1 

0.2 

1.00083 

1.20098 

1.00093 

1.20098 

2 

0.51099 

0.97138 

1.19032 

0.97138 

1.19032 

3 

0.51099 

0.97138 

1.19032 

0.97138 

1.19032 

4 

0.2 

1.00083 

1.20098 

1.00093 

1.20098 

5 

0.2 

1.00083 

1.20098 

1.00093 

1.20098 


depends on the bed wave velocity. The bed wave velocity is usually very small 
(Cunge et al. 1980) and therefore, larger A t values could be used. 

4.2.4 Aggradation Profiles In Non-rectangular Channels 

In all the previous simulations, the chanels were of rectangular shape. To demon- 
strate the applicability of the present model for non-rectangular channels, simulation 
was performed for Scenario 5. Scenario 5 considers also a five channel converging 
network as shown in Fig. 4.3. The constants B 0 , Ci, C 2 , C 3 and C 4 for the cross- 
sectional shape for channels [1],[2],[3],[4] and [5] are given in the Table 4.1 In this 
test case, Iowa model was used for computing the sediment carrying capacity. The 
initial bed slopes S 0 were taken as 0 . 0012 , 0.00082, 0.0021 for the root channels 
[ 1 ], [ 2 ] and [3], respectively. Sq value of 0.0012 was taken for the channels [4] and 
[5]. The discharge in the root channels i.e.,[l],[4] and [5] was equal to 0.0062m 2 /s. 
Lengths of the channels were equal to 30m each. Particle mean size and the porosity 
were taken as 0.32mm and 0.4, respectively. Computations were made for 3000sec 
with At = 1 . 0 sec. Figures 4.18-4.19 show the change in the bed profile in the channel 
[1] and the change in the cross-sectional shape at the inflow section after 3000sec. 
Identical results were obtained for all the root channels as expected. In the numerical 
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Figure 4.20: Temporal variation of cross section at 10m from U/S, Channel 1, Scenario- 
5 


procedure described in section 3.6.3, it was assumed that all the deposition occurs 
on the bed and bed rises in parallel layers as a function of time (Chang and Hill 
1976). 

Results presented in Fig 4.19 confirm that the above procedure has been correctly 
adopted in the computer code. There was no aggradation beyond 8m (Fig. 4.18) 
because the bed wave had not arrived at this point within the computational time 
period. It should be noted here that the bed levels and the cross-section would not 
change until the bed wave arrives from the upstream side because the initial flow 
conditions were uniform through out the network. This is also reflected in Fig. 4.20 
for the cross-section at x=10m. 
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Table 4.2: Cross-sectional shape parameters for Scenario-6 


channel no. 

Cross-sectional Shape Parameters 

Initial bed slope 
So 

B 0 (m) 

Cl 

c 2 

c 3 


2 


0.94934 

1.7714 

0.94934 

1.7714 

0.0031 

3 


0.94934 

1.7714 

0.94934 

1.7714 

0.00066 

4 


0.94934 

1.7714 

0.94934 

1.7714 

0.00066 

5 


0.94934 

1.7714 

0.94934 

1.7714 

0.00026 

6 


0.94934 

1.7714 

0.94934 

1.7714 

0.00026 


3 



4.2.5 Degradation Due To Sediment Depletion 

In the proor discussed simulations, aggradation in a converging channel system 
was considered. Applicability of the model to diverging tree-type networks and 
to simulate degradation is demonstrated in this section. In the simulation for the 
Scenario 6, a five channel diverging tree-type network as shown in Fig. 4.21 was 
considered. Length of each channel was equal to 30m. The specified values of the 
initial sectional shape parameters Bq,Cx,C 2 ,Cz,C\ and the initial bed slopes, So 
are presented in Table 4.2. Porosity of the sediment bed layer was taken as 0.4. 
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Figure 4.22: Bed profile along Channel 2, Scenario-6, Procedure 1 


The mean size of the sediment particles and the specific gravity of the particle were 
taken as 0.32mm and 2.65 respectinely. A flow depth of 0.075m was specified as 
the downstream boundary condition for all the pendant channels i.e.,[6],[5] and [3]. 
Total discharge in the root channel [2] was taken as 0.0062m 3 /s. The sediment 
carrying capacity and the friction slope were computed using the Iowa model. 
Simulation for the Scenario 6 considers the degradation of the channel system due 
to complete depletion of the sediment load into the root channel [2] (Fig. 4.21). 
Simulation were performed for a total time =3000sec. with A:r=2m, and At=1.0 
sec. While performing the simulations for the Scenario 6, two different procedures 
as discussed in chapter 3 were adopted for modelling the cross-sectional change due 
to degradation. In the procedure 1, the degradation occured without changing the 
cross-sectional shape while, in the procedure 2, the degradation was proportional to 
the local shear stress (Chang and Hill 1976) 

Figure 4.22 shows the change in the bed profiles after 3000sec. for the root 
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Figure 4.25: Temporal variation of cross section at inflow to Channel 1, Scenario- 
6, Procedure 2 


channel i.e.,[2] when procedure 1 was adopted for modelling the changes in the cross- 
section. It is obvious from this figure that the degradation wave had not arrived at 
the junction point. Consequenty, no degradation occured (not shown here) in any 
other channel of the system. Figure 4.23 shows the change in the cross-section after 
3000 seconds at the inflow section of the root channel. It can be clearly seen from 
Figure 4.23 that the procedure 1 for degradation has been correctly adopted in the 
computer code. Figures 4.24 and 4.25 show the change in the bed profile for the root 
channel and the cross-section at the inflow to the root channel after 3000 seconds 
when the procedure 2 was used to affect the change in the cross-section shape. 
A comparision of Figs. 4.22 and 4.24 indicates that the procedure adopted for 
proportioning the degradation across the cross-section affects the overall simulation 
results significantly. It can be observed that the bed wave had moved faster in the 
case of procedure 2 (Fig. 4.24). By proportioning the degraded area with respect to 
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the local shear stress, one gets more degradation at bottom. Consequently the bed 
slope for the reach changes much more than in the case of slope obtained by parallel 
degradation. Naturally these larger changes in the bed profile propagate faster. 
Once again it is emphasized here that the way to acount for change in the cross- 
section is a matter of engineering judgement. Conclusions regarding the correctness 
of a procedure should await more information from laboratory and field tests. The 
simulations presented here indicate only that the procedures are correctly encoded 
in the computer code. 
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Chapter 5 


Summary And Recommendations 
For Furthur Investigations 


This investigation presented a numerical model for the aggradation and degradation 
in tree-type channel networks. The model is based on a quasi-steady flow assumption 
and an uncoupled approach for the solution of the governing equations. The model 
can simulate the bed level variation in rectangular and non-rectangular channels. 

The numerical model was verified by comparing the computed results with, 
available experimental and other earlier numerical results. The numerical model 
simulated the bed level variation in the networks satisfactorily. However, The model 
did not give results when the Froude number was close to unity because the model 
solves the steady state energy equation for obtaining flow profiles. The model gave 
as good results as those obtained by a complete unsteady coupled model. The 
comparision between the simulated and measured aggradation in laboratory a flume 
was satisfactory. However, its applicability to natural river networks needs to be 
tested. The numerical model presented in this study is applicable for channels with 
beds composed of uniform meterial. This model dose not consider the bank stability. 
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The following is a list of subjects that require further research. 

1. Mechanism of aggradation or degradation process at the confluence points in 
networks needs to be studied in more detail. 

2. Two dimensional models need to be devoloped using quasi-steady technique. 

3. Bed material in natural channels is seldom uniform and bed degradation is 
greatly controlled by armouring. Research needs to be done for inclusion of 
the effect armouring in models for networks. 
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